if (!exists("app_root"))
    {
    stop("Need to have app_root defined to source this template");
    }

source( paste( app_root, "cgi-lib/R/ep.io.R", sep="/" ));
source( paste( app_root, "cgi-lib/R/bga_ord.R", sep="/" ));

# fix for bga - TEMP-TEMP-TEMP!!!
bga_1 = function( dataset, classvec, type = "coa", ... )
    {
    classvec<-checkfac(classvec);
    nclasses=length( levels( classvec ) );

    # Run the ordination part
    data.ord = ord( dataset, type = type, classvec = classvec, trans = TRUE );
    ord.class = class( data.ord$ord );

    ## Because between uses as.matrix, it will add X in front of any rownames
    ## starting with a numeral. Therefore tidy up
    genenam = sub( "^X", "", colnames( data.ord$ord$tab[ sapply( rownames(dataset), function(x) regexpr( "^\\d", x, perl = TRUE ) == 1 ) ] ) );
    if (!is.null(dim(genenam)))
        {
        colnames(data.ord$ord$tab) = rownames(data.ord$ord$co)= rownames(data.ord$ord$c1) = genenam;
        }

    # Run Between Group analysis, and return class dudi.bga
    data.bet = between( data.ord$ord, classvec, scannf = FALSE, nf = nclasses - 1 );

    res = list( ord = data.ord, bet = data.bet, fac = classvec );
    class(res) = c( type, "bga" );
    return(res);
    }


ep.run_bga = function( src, testsrc, filename, testfilename, outputdir, classvec, ordtype, outputlist, supplproj, supploutput )
    {
    outbase=paste(outputdir,"/",filename,sep="");

    ################################
    # Read data
    ################################
    dataset = ep.readBin(src);
    classvec = strsplit(classvec, ",")[[1]];

#   colnames(dataset)<-srccolannot

    ################################
    # Run BGA
    ################################
    otype= tolower(ordtype);

#   if ( !supplproj=="Y" )
        bet.res = bga_1( dataset, classvec, type = otype );




    ################################
    # Project suppl data - TO DO
    ################################

#   Run BGA+ Supplementary Analysis!
#   if ( supplproj == "Y" )
#       {
#       print(paste("Running Suppl Analysis on", testfilename, sep=" "));   
#       testdata <- ep.readBin ( testsrc );
#       bet.res<-bga.suppl(dataset, supdata=testdata, classvec=classvec);
#       sup.res<-suppl(bet.res, testdata);
#       outputlist = paste(outputlist,supploutput, sep="");   
#       }

    ################################
    # Produce output summary
    ################################
    sink( paste( outbase, "_log.txt", sep = "" ) );
    print( bet.res );
    sink();

    ################################
    # BGA output files
    ################################
 
    outputlist = strsplit( outputlist, "" )[[1]];

    sink( paste( outbase, "_list.txt", sep = "" ) );
    print(outputlist);
    sink();

    ################################
    # Write file, eigenvalues,
    # column (array),
    # and row (lines, genes) coordinates
    ################################

    write.res = function ( res = bet.res$bet, fileExt = "_eig.txt", outdata = "eig" )
        {
        filenam = paste( outbase, fileExt, sep = "" );
        data = signif( res[[outdata]], 3 );

        # if data has more than 300 rows only retain first 100 and last 100
        # TO DO - disabled?
        if ( 1 == 0 && is.data.frame(data) && nrow(data) > 300 )
            {
            p = nrow(data) - 100;
            data = data[-101:-p,];
            }

#       if ( nrow ( data ) > 300 ) \{
#           p <- nrow(data) - 100;
#           data <- data[-101:-p,]
#       \}

        write.table( data, sep = "\t", quote = FALSE, row.names = TRUE, col.names = FALSE, file = filenam );
        }

    write.res( outdata = "eig", fileExt = "_eig.txt" );

    ################################
    # Define plot
    ################################

    plot.res =function( ord, fileExt = "_array.png", type = "A" )
        {
        filenam = paste( outbase, fileExt, sep = "" );
        png( filenam, width = 600, height = 600 );
        if ( type == "D") plot(ord);
        if ( type == "A") plotarrays(ord);
        if ( type == "G") plotgenes(ord);
        if ( type == "E") scatterutil.eigen(ord$eig); 

        if ( type == "S") plotarrays(ord);

        dev.off();
        }

#   plot.res( ord = bet.res, fileExt = "_array.png", type = "A" );

    write.pdb = function( ord, fileExt = "_li", coord = "li", classv = NULL )
        {
        filenam = paste( outbase, fileExt, sep = "" );
        xyz = ord[[coord]];
        if ( ncol(xyz) > 2 )
            html3D( xyz, classvec = classv, writepdb = TRUE, scaled = TRUE, filenamebase = filenam );
        }

    write.pdb( bet.res$bet, coord="ls", classv = bet.res$fac );

    ################################
    # Parse output wish list
    ################################

    for ( i in outputlist )
        switch( EXPR = i,
            # Output files default plot (D), eigenvalues $eig (E), array $ls (R),
            # gene $co (C) and centroid $li (T) co-ordinates
            E = write.res( outdata = "eig", fileExt = "_eig.txt" ),
            R = write.res( outdata = "ls", fileExt = "_ls.txt"  ),
            C = write.res( outdata = "co", fileExt = "_co.txt" ),
            T = write.res( outdata = "li", fileExt = "_li.txt" ),
            P = write.table( bet.res$suppl, file = paste( outbase, "_suppl.txt", sep = "" ), sep= "\t" ),
            # Default plot
            D = plot.res( ord = bet.res, fileExt = "_def.png", type = "D" ),
            # Graph of Eigenvalues
            V = plot.res( ord = bet.res$bet, fileExt = "_eig.png", type = "E" ),
            # Graph of Arrays
            A = plot.res( ord = bet.res, fileExt = "_array.png", type = "A" ),
            # Plot of genes
            G = plot.res( ord = bet.res, fileExt = "_genes.png", type = "G" ),
            # Plot of supplementary data 
            S = plot.res( ord = bet.res$suppl[,1:2], fileExt = "_suppl.png", type = "S" )
            );


        print( paste( "Graphs and output files in ", outputdir, "/", filename, sep = "" ) );
    }
